% Plot the power at M2 (12.421 hours) of the observatory data on to the
% predicted power map

colormax = 10;
%S = dir('/home/mnair/projects/tides/*spectra_new.mat'); % Only PMTM results (variable P) is useful
S = dir('/home/mnair/projects/tides/*_1995_2010_ROB_DAY_NIGHT_7_years_SWEEP.mat');
for i = 1:length(S),
    robstn(i,:) = S(i).name(1:3);
end;
%load '/home/mnair/projects/tides/M2_Abs_Colormap' map1;
%open('/home/mnair/projects/tides/M2_Bz_Abs_0_720_1440.fig');
open('/home/mnair/projects/tides/M2_Hr_power_720_1440.fig');
set(gca,'FontSize',16);
colorbar('h');
load /home/mnair/projects/tides/M2_results_720_1440 Hr_power;
caxis([0,colormax])
map1=colormap;
hold on
periods = [4 4.8 6 8 11.967236 12  12.421 12.6583 23.934472 24 25.891];
% period_string = ['S6';'S5';'S4';'S3';'K2';'S2';'M2';'N2';'K1';'S1';'O1'];

%OK I don't need periods < 10

periods_12 = [11.967236 12  12.421 12.6583 ];
periods_12_text = [11.9 12.01  12.43 12.67 ];
period_string_12 = ['K2';'S2';'M2';'N2'];

periods_24 = [23.934472 24 25.891];
periods_24_text = [23.8 24.01 25.9];
period_string_24 = ['K1';'S1';'O1'];

% 
% diff_ls_all = 0;
% diff_ls_nigh = 0;
% diff_rob_all = 0;
% diff_rob_night = 0;

nstsn_plot = 0;

%fid = fopen('load /home/mnair/projects/tides/M2_results_table.txt','wt');

m2_global_power = zeros([1,100]);
for nstn = 1:length (S);
    
    eval(['load /home/mnair/projects/tides/' S(nstn).name ]);
    
    
    
%    median_amp = median(sqrt(P/(en-st)));
 median_amp = 1;
    
    if median_amp < 1.5
        
            
%            L = 1./(F*3600) >= 12.419 & 1./(F*3600) <= 12.4240;
            
%             if any(L)
%                 m2_power_pmtm = mean((P(L,1)/(en-st))); %This power doesn't mean anything
%             end;
%             
            % get the m2 power from different processing methods
%             m2_power_night   = spectra_data_night_ls(7)^2;
%             m2_power_all = spectra_data_all_ls(7)^2;
%             m2_power_all_rob = spectra_data_all_rob(7)^2;
%             m2_power_night_rob = spectra_data_night_rob(7)^2;
            
                        % get the m2 power from different processing methods
            m2_power_night   = spectra_data_day_rob(7)^2;
%             m2_power_all = spectra_data_(7)^2;
%             m2_power_all_rob = spectra_data_all_rob(7)^2;
            m2_power_night_rob = spectra_data_night_rob(7)^2;

            
            % Average the predicted power over a square
            st_lat = round((90-latitude)*4)-2;
            en_lat = round((90-latitude)*4)+2;
            st_lon = round(longitude*4)-2;
            en_lon = round(longitude*4)+2;
            
            if st_lat < 1,
                st_lat = 1;
            end;
            
            if en_lon > 1440,
                en_lon = 1440;
            end;
            
             if en_lat > 720,
                st_lat = 720;
            end;
            
            if st_lon < 1,
                st_lon = 1;
            end;
            
            pred_m2_power = mean(mean(Hr_power(st_lat:en_lat,st_lon:en_lon))); %Average over an area
            pred_m2_power_1cell = Hr_power(round((90-latitude)*4),round(longitude*4)); % Value in the cell
            
            
            % calculate  error
            nstsn_plot = nstsn_plot + 1;
%             diff_ls_all(nstsn_plot) = pred_m2_power - m2_power_all ;
%             diff_ls_night(nstsn_plot) = pred_m2_power - m2_power_night ;
%             diff_rob_all(nstsn_plot) = pred_m2_power - m2_power_all_rob ;
%             diff_rob_night(nstsn_plot) = pred_m2_power - m2_power_night_rob ;
%             
%             diff_ls_all_1cell(nstsn_plot) = pred_m2_power_1cell - m2_power_all ;
%             diff_ls_night_1cell(nstsn_plot) = pred_m2_power_1cell - m2_power_night ;
%             diff_rob_all_1cell(nstsn_plot) = pred_m2_power_1cell - m2_power_all_rob ;
%             diff_rob_night_1cell(nstsn_plot) = pred_m2_power_1cell - m2_power_night_rob ;
% 
            % relative difference
            
%             rel_diff(nstsn_plot) = abs( (diff_rob_night(nstsn_plot) - pred_m2_power ) / pred_m2_power );
            % Find color to plot
            
            %cidx = floor(63*(m2_power_night)/(colormax)+1); % M2 LS NIGHT
            cidx = floor(63*(m2_power_night_rob-m2_power_night)/(1)+1);
            
            if cidx > 64,
                cidx = 64;
            elseif cidx < 1,
                cidx = 1;
            end;
            
            figure(1);
            plot(longitude, 90-latitude, 'o', 'MarkerFaceColor',map1(cidx,:),'MarkerSize',20, 'MarkerEdgeColor','white');
            %end;
            
            
%             figure(2);
%             
%             set(gcf,'Position',[ 680   382   887   713]);
% 
%             
%             subplot(211);
%             
%              T = (1/60) * (en-st) ; % band width of FFT
%             
%             L = 1./(F*3600) >= 11.5 & 1./(F*3600) <= 13;
%             plot(1./(F(L)*3600), log10((P(L)/T)), 'b', 'LineWidth', 2);
%             set(gca,'FontSize',16);
%             hold on;
%             text(periods_12_text,zeros([1,length(periods_12)]),period_string_12);
% %             plot(periods(7),log10(pred_m2_power),'rs','MarkerSize',10);
%            
%                 %plot(periods(5:8), log10(spectra_data_all_ls(5:8).^2),'r*');
% %                 plot(periods(5:8), log10(spectra_data_night_ls(5:8).^2),'k*');
% %                 plot(periods([5,7,8]), log10(spectra_data_night_rob([5,7,8]).^2),'ko');
%          
%             
%             legend('pmtm','MODEL','LS NIGHT','RB NIGHT');
%             a = axis;
%             
%             for kk = 1 : length(periods_12),
%                 line([periods_12(kk) periods_12(kk)], [a(3) a(4)],'LineStyle','--', 'color','k');
%             end;
%             
%             
%             title(['Spectra at ' S(nstn).name(1:3) sprintf(' (lat=%5.2f,lon=%5.2f )'...
%                 ,latitude,longitude)]);
%             xlabel('Period - Hours');
%             ylabel('log(nT^2)');
%             hold off;
%             
%             subplot(212);
%             L = 1./(F*3600) >= 23.5 & 1./(F*3600) <= 26.5;
%             plot(1./(F(L)*3600), log10((P(L)/(en-st))), 'b', 'LineWidth', 2);
%             set(gca,'FontSize',16);            
%             hold on;
%             text(periods_24_text,zeros([1,length(periods_24)]),period_string_24);
%             
%        
%                 %plot(periods(9:11), log10(spectra_data_all_ls(9:11).^2),'r*');
% %                 plot(periods(9:11), log10(spectra_data_night_ls(9:11).^2),'k*');
% %                 plot(periods(9:11), log10(spectra_data_night_rob(9:11).^2),'ko');
%        
%             title(sprintf('( pred m2=%5.2f, night m2 ls=%5.2f, night m2 rob =%5.2f )',...
%                 pred_m2_power, m2_power_night,m2_power_night_rob ));
%             a = axis;
%              h = legend('pmtm','LS NIGHT','RB NIGHT');
%             for kk = 1 : length(periods_24),
%                 line([periods_24(kk) periods_24(kk)], [a(3) a(4)],'LineStyle','--', 'color','k');
%             end;
%             
%             xlabel('Period - Hours');
%             ylabel('log(nT^2)');
%             %%
%             eval(['print -dpng  /home/mnair/projects/tides/' S(nstn).name 'power_spec' ]);
%             hold off;
%             pause;

% just to see how the global power of M2 look like

m2_global_power = m2_global_power + spectra_data_rob_m2_n;
        end;
    end;
